Topological order and the deconfinement transition in the (2+1) dimensional compact 
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We study an Abelian compact gauge theory minimally coupled to bosonic matter with charge q, 
which may undergo a confinement -deconfinement transition in (2+l)D. The transition is analyzed 
using a nonlocal order parameter W, which is related to large Wilson loops for fractional charges. 
We map the model to a dual representation with no gauge field but only a global q-state clock 
symmetry and show that W correspond to the domain wall energy of that model. W is also directly 
connected to the concept of topological order. We exploit these facts in Monte Carlo simulations to 
study the detailed nature of the deconfinement transition. 
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I. INTRODUCTION 



The appearance of excitations with fractional quan- 
tum numbers in strongly correlated systems has received 
much interest in recent y ears 4iSiMi5iSiLSi2iifi The splin- 
tering of the electron into a neutral spinon and a charged 
spinless holon, and the resulting spin-charge separation, 
was proposed early after the discovery of the supercon- 
ducting high-Tc cuprates as a route to explain the non- 
Fermi liquid behavior observed in these systems. The 
fractionally charged Laughlin quasi-particles in the frac- 
tional quantum hall effect is another famous example. 
While spin-charge separation is well established in one 
dimensional electron gases, the occurrence in higher di- 
mensions remains poorly understood. Usually, the frac- 
tionalization in dimensions larger than one, is imag- 
ined to occur via a confinement-deconfincmcnt transition 
of an effective low energy gauge theory description of 
the system^i*i£ A number of such gauge theories have 
been proposed for the underdoped cuprates, including 
E/(l),W, 4 ,Wi 8 i 13 , 14 SU(2)jL and Z^gauge theories. 
In addition, the possibility of fractionalized phases in 
models containing only bosons has been proposed^ The 
fractionalized phases are not ordered in a conventional 
sense, which breaks any obvious symmetries of the sys- 
tem. Instead they are characterized by a special kind 
of order called topological order>i& Topological order is 
robust against local perturbations and furthermore leads 
to a ground state degeneracy. It has been proposed that 
this topologically protected degeneracy could serve as a 
building block for qubits in a quantum computer An 
important prediction is that in finite sized systems the 
degeneracy is lifted due to tunneling processes, leading 
to a splitting A ~ e" i/? J£ 



In this paper we study the relation between topologi- 
cal order and deconfinement in one of the simplest mod- 
els displaying topological order — the compact Abelian 
Higgs model. It should be noted that a realistic the- 
ory of the high-T c cuprates would require also the in- 
clusion of gapless nodal quasiparticles, which are here 
neglectedi 19 i 20 i 21 i 22 i 23 We thus consider a (2+l)D com- 



pact U(l) gauge theory minimally coupled to a charge-q 
bosonic matter field, with an Euclidean action 

S = - Jj^ cos(V M r - q A r ^) C0S ( fi r J, (1) 



r// 



where 9 r and A rpi are compact phases € [0, 27r) liv- 
ing on the sites and links of a 3D simple cubic lat- 

efj, u \V v A r x is the dual field 



tice, respectively. B 
strength, with the lattice difference operators defined by 
V M / r = V M / r+ e (l = / r +e^ - /r- The pure gauge the- 
ory in absence of matter is always confining in (2+l)D, 
due to the proliferation of instantonsi 2 ^ (In (3+l)D a 
deconfined Coulomb phase is also possible.) 2 ^ Partly be- 
cause of this much attention has focused on Z2-gauge 
theories^ which do allow a confinement-dcconfincment 
transition in (2+l)D. Alternatively, when the U (1) model 
is coupled to a matter field which does not belong to the 
fundamental representation (i.e., when the charge q of 
the boson differs from unity in Eq. l|T|ll. a phase transi- 
tion is possible^ This could happen, e.g., if the bosonic 
field describes the pairing of spinonsASSi Indeed, the U(l) 
gauge theory coupled to a charge 2 boson is equivalent 
to a Z 2 -gauge theory in the limit J — > oo. In such a 
system, the fractionalized phase would then correspond 
to a phase where the spinon-pairs breaks up into sep- 
arate free excitations. More generally, the model with 
gauge charge q reduces to a Z g -gauge theory in the limit 
J — * oo. In the opposite limit g — > 0, the gauge field fluc- 
tuations freeze out and a 3D XY model (with a global 
U(l) symmetry) results. One may ask how the transition 
interpolates between these two extremes with increasing 
screening length A = 1/ yj Jgq 2 . The naive expectation 
would be that the A — > oo transition would represent an 
unstable fixed point and that any finite value of A would 
flow to zero upon renormalization, giving a transition in 
the Z q universality class. In this case the 3D XY point at 
A — > oo would be an isolated point in the phase diagram. 

Recently, large scale Monte Carlo simulations of the 
(2+l)D compact U(l) Higgs model were performed by 
Sudb0 et alwL and Smiseth et al£& for several values of 
the charge q > 1 . The results of these indicated that the 
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phase diagram might be surprisingly complicated, with 
continuously varying critical exponents along the transi- 
tion line between the Z q and 3D XY values. This would 
imply that the transition line for intermediate A might 
be described by a line of fixed points. These interesting 
results were obtained from a finite size scaling analysis of 
the third moment of the action, which is a purely ther- 
modynamic observable. Although their method seems 
to work very well, it would be desirable to have a clear 
order parameter which distinguishes the confining and 
deconfining phases, especially in light of their surprising 
results. 

In this paper we employ a nonlocal order parameter W, 
previously introduced in Refl2^. which provides a direct 
probe of the confining properties of the theory. We show 
that this order parameter can be mapped via a duality 
transformation to the domain wall energy of a Z q clock 
model, and present an in depth discussion of its relation 
to top olog ical order, including many of the details left out 
in Refl2iA We obtain the finite size scaling properties of 
W, discuss its relation to the splitting A of the nearly 
degenerate ground states, and the energy gap to vortex 
excitations in the deconfincd phase. Finally, we use it in 
Monte Carlo simulations to study the detailed nature of 
the deconfinement transition of the model, and compare 
with the earlier results mentioned above. 

The paper is organized as follows. In Sec. Ill Al we re- 
formulate the theory in terms of its topological defects, 
which arc vortices and instantons, and discuss how W can 
distinguish the different phases. Going over to another 
dual formulation a model with a global g-state clock sym- 
metry (and no gauge fields) obtains, and W translates 
into the domain wall energy of that model, or equiva- 
lently the energy cost for twisting the boundary condi- 
tions by 27r/g, see Sec. Ill Bill CI In Sec. IIIDl we discuss 
the relation to topological order, addin g de tails and gen- 
eralizations of the ideas discussed in Refl29l In Sec. lIIII we 
describe Monte Carlo simulations of the (2+l)D compact 
Higgs model carried out in the loop-gas representation of 
flux lines and monopoles. Our results and conclusions 
are discussed in Sec. II VI andlVl 



II. DUALITY TRANSFORMATIONS, WILSON 
LOOPS AND TOPOLOGICAL ORDER 

In order to transform Eq. Q to alternative representa- 
tions it is convenient to instead study the Villain version 
of the model, which allow the action to be transformed 
to dual representations exactly, without further approx- 
imations, while retaining the important physics. In the 
Villain approximation the action is 




(2) 

where n rfl and k rfJi are dummy integers to be summed 
over. The summation over k ru has the effect of mak- 



ing the energy cost of putting an integer multiple of q 
flux quanta $o = 27r/q through a plaquette of the lat- 
tice vanish. The net effect of this is the proliferation of 
Dirac strings, which leads to the appearance of magnetic 
monopoles of charge q$o (= 27r). 

A standard set of transformations'^*^ allows Eq. 
to be rewritten in terms of flux lines and monopoles, in- 
teracting via the action 

K KX 2 a 2 

s = J2Y TaT ' Vrr ' inr ' + ~^~ NrVTr ' Nl '- (3) 

rr' 

Here m r 6 Z 3 is the vorticity on the links and iV r e Z is 
the monopole charge^ and K = (2ir) 2 J. The flux lines 
are constrained to start or end on the monopoles only 
in quanta of q, i.e., V • m r = qN r . The flux lines and 
monopoles interact through a screened coulomb interac- 
tion, defined by its Fourier transform V^ -1 = k 2 , + A~ 2 , 
where n 2 = J2v K kv an d K k</ = 2sin(/c„/2), i.e., 

* = i?i: (4) 

k K 

which for large distances is a screened coulomb interac- 
tion V(r) = e~ r l x /4iTr with a screening length A given 
by A~ 2 = Kg/® 2 , = Jgq 2 . In the limit A — > the in- 
teraction reduces to an onsite interaction V r = A 2 <5 r , and 
the action becomes 




The constraint V • m = qN means that a +(— ) 
monopole has precisely q outgoing (incoming) flux lines. 
For q > I one may then envisage two distinct phases 
of the system. Either the monopoles mostly pair up in 
neutral pairs, where the q outgoing flux lines of the plus 
charge all end at the minus charge of the pair, as illus- 
trated in Fig. |TTa). This phase will be realized in the 
limit of large K, where the flux lines cost a lot of energy. 
As discussed in the next section, this magnetically con- 
fining phase corresponds to deconfined fractional electric 
charges. The other possibility is that the outgoing flux 
lines mostly end on different monopoles, forming a large 
connected tangle of flux lines and monopoles, which per- 
colate through the whole system, as in Fig.QJb). Recent 
simulations' 31 support this interpretation. A transition 
between these two geometrically different phases is ex- 
pected at some critical value of the coupling K c (or g c ). 
This picture can be made more precise by studying the 
Wilson loop for fractional charges. 

A. Wilson loops 

The presence of percolating flux lines in the system can 
be detected by counting the number of flux lines crossing 
a given surface S that form a cross-section of the system, 
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FIG. 1: (Color online) Possible phases of the system. To the left the topological ordered, deconfined phase, where vortices 
cost much energy. To the right the percolating vortex tangle in the confined phase. The large Wilson loop W surrounds the 
shaded area as indicated in the figures. 



see Fig.^ We assume here periodic boundary conditions 
for the vorticity m in all directions, i.e., the system has 
the topology of a 3-torus. This corresponds to allowing 
twists in the boundary conditions for A in Eq. or 
@p*ii Let the winding number M v be the number of such 
flux line crossings (counted with signs) for a cross-section 
perpendicular to the indirection, i.e., 

M v = / m • dS. (6) 
Js 

(For notational convenience we use a continuum formu- 
lation here.) The precise location of the surface does not 
matter. In the paired phase this number will be an inte- 
ger multiple of q, whereas in the percolating phase it can 
be any integer. It turns out that precisely this property is 
captured by a large Wilson loop for fractionally charged 
test particles, 

W = W(C) = (exp (i j> A • dr\ \ , (7) 

where the loop C encircles the surface S. The expression 
in the exponent, 

/ A • dr = [ B • dS = — [ m • dS, (8) 
Jc Js 9 Js 

equals 2n/q times the winding number M u , and the Wil- 
son loop for a loop with this special geometry is thus 
given by 

W=(exp(^M v }\. (9) 



This becomes one in the paired phase where M v is an 
integer multiple of q [Fig.^a)]. In the percolating phase 
[Fig. mb)] no particular value of the winding number 
M v will be favored, making the W approach zero with 
increasing system size. Note that it is the global topo- 
logical structure of the flux line tangle that determines 
the value of W, while small scale entanglement is unim- 
portant. 

A Wilson loop decaying with an area law, W(C) ~ 
exp(— cL 2 ) implies a linearly confining attraction be- 
tween static fractional charges, while a perimeter law, 
W(C) ~ exp(— cL), means that such charges would be de- 
confined. Note that the behavior of the charges is oppo- 
site of the monopoles: The charges are deconfined when 
the monopoles form pairs and confined when the pairs 
break up. Remember that the charge of the bosons in 
the original model is q, so the Wilson loop for charge-g 
test particles would include a factor q in the exponent, 
making W equal to 1 exactly in both phases. For frac- 
tional charges, however, it changes from unity to an area 
law W ~ exp(— cL 2 ), when going through the transition 
(deep in the confined phase one finds c = g/2). Another 
way to say this is that the charge-g dynamic matter can- 
not screen out the fractionally charged test particles. 

The absence of a perimeter law for W in the decon- 
fined phase is a consequence of the special geometry of 
a loop which covers a whole cross-section of the system, 
and in some sense has no perimeter due to the periodic 
boundary conditions. This makes our order parameter 
W different from ordinary Wilson loops, and is a big 
advantage since it makes it easier to analyze its scaling 
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behavior as discussed in Sec. Ill CI For a finite sized loop 
in a much larger system the Wilson loop would follow a 
perimeter law rather than approach unity in the decon- 
fined phase. This is due to configurations where part of 
the flux connecting a monopole-anti-monopole pair would 
pass through the loop and part outside the loop. In the 
paired phase this can happen only very close to the loop 
itself, hence giving rise to a perimeter law. In the next 
two sections we describe the relation of the Wilson loop 
order parameter W to the domain wall energy in a dual 
clock-symmetric model. 



B. Mapping to a model with clock symmetry 

We will now map Eq. J3J to a dual model with a g-state 
clock symmetry. The idea is to start from the action 

K K\ 2 n 2 
S = ^2—m T -V TT >nL r , + —?-N r V Tr >N r , 



2 2 

rr' 

i-Xr (V • m r - qN r ) . 



(10) 



where the constraint V • m = qN is implemented via a 
Lagrange multiplier Xr, and then integrate out the flux 
lines and monopoles to get an effective action for x- We 
decouple the vortex-vortex and monopole-monopole in- 
teractions with the help of two auxiliary fields, a and <j), 
to get 

rr' 

+ ia r ■ m r + i(j) r qN r + i\ r (V • m r — qN T \ll) 

r 

Now the summation over the integers m r and N r can be 
performed, leading to 

a = Vx + 27rn, ^ — x + ^nk/q, n, k integers. (12) 

We then obtain the action 

s = tk (VXr + 27rrir) v * {VXr ~ + 27air,) 



2KX- 



( X r + 2wkr/q) V-} (Xr- + 2nk r , / q) . (13) 



As announced this model has only a global g-state clock 
symmetry and contains no gauge fields and only short 
range interactions. Rewriting it in Fourier space gives 

S = ^E( K2 + A_2 )l^ + 27rn k | 2 



2KQ 

k 

1 + A" 2 ) |Xk + 27Tk k /q\ 2 . 



(14) 



Obviously 1 = n 2 + A 2 is short ranged and may at 
least for small values of A (i.e. large g) be replaced by 



a completely local interaction, V^} « \~ 2 8 rr i. The cor- 
rections enter with higher order derivative and should 
therefore be irrelevant at small A. Dropping these and 
"undoing" the Villain approximation now gives 



I r/i r ) 



(15) 



which clearly has a g-state clock symmetry. For small A 2 
it is legitimate to replace the continuous variable Xr by 
a discrete one Xr ~ 2irk r /q leading to 



1 



A'A 2 



cos(27rV,iA: r /g), 



(16) 



which is just the familiar g-state clock model. Equa- 
tion (|16(l is exact in the limit A — >• 0, where the well 
known duality between Z q -gauge theory and the g-state 
clock model is regained. Equation (|13|l remains valid for 
arbitrary A. 

Note that for fundamental matter, when q = 1, 
Eq. (|13[) essentially describes an XY model in a finite 
external magnetic field, which does not have a phase tran- 
sition; this agrees with Ref. 0. 



C. Domain wall energy 

We now proceed to show that the domain wall en- 
ergy in the g-state clock symmetric model, Eq. <|13ll . 
corresponds precisely to the large Wilson loop W in 
the original compact gauge theory. Consider a finite 
sized system with side length L and periodic bound- 
ary conditions. A domain wall can be forced into the 
system by changing the boundary conditions such that 
x(r + Le, y ) = x(r) + Ax, Ax = 2irk/q. The domain wall 
energy is defined as the free energy cost AF for introduc- 
ing such a twist. It is convenient to make a change of vari- 
ables in the twisted system to get back to periodic bound- 
ary conditions by setting x( r ) — $( r ) + A0(r), where 9 
obeys periodic boundary conditions, and A9(r) = con- 
stant except across a surface S occupying a cross section 
of the system perpendicular to the direction where 
A0(r) jumps discontinuously by 2nk/q. We can then 
evaluate the domain wall energy as 



AF = F t 



twisted 



F 



untwisted 



= -ln< 



e AS ) 



(17) 



where the average (• ■ • ) is taken with respect to the un- 
perturbed system, and AS = S[9 + A9] - S[9}. Going 
back to the flux line representation Eq. (|10|) we have 

AS = ^A6» r (V-m r -gA r ) 

r 

= -i ^2 m r • VA6» r + A9 r qN r 



2irik 



Y m ' n r - i Y &9 r qN r (18) 



r£5 



where S is the surface at which Ad is discontinuous and 
n r is the surface normal at lattice site r. Because A8 r 
is an integer multiple of 2n/q only the first term on the 
right hand side contribute when exponentiated. The re- 
sulting formula is more transparent when written in the 
continuum, 

( e - AS ) = (e*?* = (e^ M ») . (19) 

This is nothing but the large Wilson loop W for a test 
particle with fractional charge k. Thus the domain wall 
energy of the dual model is related to W by 

AF = -\nW. (20) 

In the (dual) ordered phase this will be proportional to 
the surface area, AF ~ L^ 1 = L 2 , consistent with an 
area law for the Wilson loop W ~ exp(— cL 2 ). In the 
disordered phase, creating a domain wall costs no energy 
and AF — > with system size. Precisely at a critical 
point AF — > a universal constant, independent of system 
size. For a finite sized system we expect the finite size 
scaling relation 

AF{5,L) = f{8L l '») (21) 

to hold, where 8 = (K^ 1 — K^ 1 )/ F" 1 is the tuning 
parameter of the transition and v the correlation length 
exponent. 

From simulation data the critical point can thus be 
located by plotting AF (or equivalently W) vs the tun- 
ing parameter for different system sizes and observing 
the point where the different curves cross. We will show 
examples of this below in Fig. [3] 

D. Topological order 

One characteristic feature of systems displaying topo- 
logical order is a ground state degeneracy which depends 
on the topology of the space on which the model is de- 
fined. A further prediction is that for finite sized sys- 
tems the exact degeneracy will be lifted due to tun- 
neling events, leading to a ground state splitting A ~ 
exp (— £■/£). 16 In Ref. |2^ it was shown that there is a 
close relationship between the lifting of the ground state 
degeneracy and the large Wilson loop W discussed above. 
Here we give more details on this relation and discuss the 
generalization to topologies other than the torus. 

To simplify the discussion we consider first the case 
q = 2 on a torus, with W in the space-time plane. In 
order to study ground state properties of the system we 
have to consider anisotropic systems with size L x L x (3 
and take the limit where the length of the time dimension 
(3 goes to infinity. In the infinite system size limit the de- 
generate ground states may be characterized by having 
an even or odd number of flux quanta $0 through the 
holes of the torus, giving a degeneracy of 2 2 = 4 since 
a torus has two holes^ The general case of arbitrary q 



5 



FIG. 2: A flux quantum tunnels out from the hole of a torus. 

on a manifold with genus g (i.e., which has g handles 
and 2g holes) naturally leads to q 2g degenerate ground 
states. For finite system size the degeneracy will be lifted 
by vortex tunneling processes, in which a vortex tunnels 
around a nontrivial path on the torus as illustrated in 
Fig. El In the deconfined phase such processes are ex- 
ponentially suppressed and the system is topologically 
ordered. In the confinement phase, however, the vortices 
condense, and the tunneling is no longer negligible even 
in the thermodynamic limit. Let us concentrate on just 
one of the two holes of a torus for now. We will de- 
note the quantum mechanical states with an even or odd 
number of flux quanta through the hole by |0) and |1), 
respectively. 

The space-time configurations of vortex lines and 
monopoles that contribute to the partition function splits 
into different topological sectors, depending on whether 
the winding number M v defined in JSJ is even or odd, 
i.e., Z = Z evea + Z odd . Our Wilson loop W (for q = 2) 
can be written in terms of these as 

W = (e i7rM ») - Zcvcn ~ Zodd . (22) 

Quantum mechanically the configurations in the even sec- 
tor contribute to the diagonal matrix elements 

Z ovcn = (0| e -^|0) = (l|e-^|l), (23) 

while the odd configurations contribute to tunneling ma- 
trix elements 

^odd = (0| e -^|l) = (l| e - /3d |0). (24) 

The tunneling lifts the degeneracy and mixes the different 
flux states. The Hamiltonian restricted to the (almost) 
degenerate ground state Hilbert space is 

*-(**) < 25 » 

which has energy eigenvalues E = Fo=fA and eigenstates 
|±> = (|0) ± |1)) /\?2. The ground state splitting A can 
be related to W by the following equality, valid for j3 — > 

00: 

e -/3H _ e -f3E ( cosh /3 A sinh/3A \_( Z cvcn Z odd \ 
^ sinh/3A cosh/3A J y Z odd Z cvcn ) 

(26) 

Using Eq. (|2*2*|) we have 

W = e" 2/3A , (27) 
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as found in Ref. 29. In the deconfined phase W = 1 in the 
thermodynamic limit and we must have A = for [3 — > 
oo, i.e., an exact degeneracy. For finite L the splitting can 
be evaluated directly deep in the deconfincmcnt phase, 
where the flux lines cost much energy and are very dilute. 
The dominating configurations which have to be taken 
into account then are those with n straight tunneling 
trajectories across the system with action S ~ aLn. This 
gives 



Z — ✓ n I 



— aLn pftLe 



and 



WZ a Y, ^—j— {(3L) n e-° Ln = e~^ Le ~ 



(28) 



(29) 



where a is the line tension and there are (3L places to 
put the tunneling trajectory. In this limit W will thus be 
given by W w exp (— 2/3Le~ CTi ) . Precisely such an ex- 
ponential form of the ground state splitting has been ar- 
gued for by Wen*S Approaching the phase transition the 
flux line will no longer be straight and interaction effects 
between different flux lines must be taken into account, 
making it much more difficult to estimate the ground 
state splitting. However, matching to the scaling form 
Eq. (JUJ) for W suggests that W = exp (~2(3L/^ 2 e- aL ^) , 
so that 



L/e 



(30) 



where £ ~ \$\~ v is the correlation length. This should 
hold for L > £ on the deconfined side of the transition. 
Indeed, we will verify in the next sections that this form is 
in agreement with the data from numerical simulations. 
In the confinement phase an area law W <~ e~ c/3L for the 
Wilson loop translates into a splitting A ~ L that grows 
with system size. 

The arguments given above naturally generalize to 
manifolds of arbitrary genus g and charge q. There 
are 2g nontrivial loops around which vortex tunneling 
can take place. Each such tunneling event changes the 
flux through the corresponding hole by one flux quan- 
tum $o. In the infinite system size limit the q 2g flux 

states = |/ 2 ) . . . |/2g), where G Z q is the num- 
ber of flux quanta through hole i (mod q), are degenerate 
ground states. The tunneling lifts the degeneracy leading 
to a splitting 



(31) 



from the ground state. Here we generalize Eq. 10 to 



^=^exp(^X>M> 



(32) 



where Mj measures the flux quanta (mod q) tunneling 
through the space-time surface surrounding hole i, and 



ki G Z q . The corresponding eigenstates are given in the 
appendix in Eq. (|A.2|I and the ground state is the product 
of the symmetric combinations of all the flux states 



\G) 



29 / q-l 

n MfEiA 



(33) 



Away from the transition one may, to a very good approx- 
imation, neglect higher order tunneling processes and 
only consider those which tunnel single flux quanta. De- 
noting those tunneling matrix elements by Tj the energy 
levels become E? = E — ^ i T i 2cos{2TTk i /q), i.e., the 
contributions from individual tunneling processes are just 
added together. 



E. Vortex gap 

An important property of the ordered phase is that 
vortex excitations are gapped. This gap can be related 
to the large Wilson loop oriented normal to the time 
direction W T . Returning to the q = 2 case the vortex 
gap is 

E v =- hm i In = - l im I m Ll^L , (34) 

P^c* (3 Z cven P^oo (3 1 + W T 

where Z dd and Z even now refer to the partition functions 
with an odd and even number of vortices present. A 
scaling argument analogous to Eq. (p?TTf> gives 



Sin 4, 



(35) 



on the deconfined side of the transition in the low tem- 
perature limit — > 00 and £ < L. 

III. MONTE CARLO ALGORITHM 

We use Monte Carlo simulations to obtain quantitative 
results for the model and the nonlocal order parameter 
W. The simulations are performed in the loop-gas rep- 
resentation, Eq. Q, of flux lines and monopoles, and we 
focus on the cases q = 2 and q — 3. We use simple cubic 
lattices with periodic boundary conditions for the vor- 
tices and monopoles, corresponding to fluctuating twist 
boundary conditions in the phase representation. 29 

In the simulation we must update both the flux lines 
and the monopoles. In order to calculate W it is nec- 
essary to include global moves that change the winding 
numbers M v and thus the topology of the vortex con- 
figuration. In a conventional Metropolis algorithm the 
acceptance ratio for such moves becomes exponentially 
small with increased lattice size. To overcome this we 
use a worm cluster Monte Carlo algorithmic^ This al- 
gorithm is described for on-site interactions in Ref. I34L 
but is straightforward to adapt to longer range interac- 
tions. 
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FIG. 3: The large Wilson loop order parameter W as function of coupling constant, for different values of the screening length 
A. Upper left: W vs g/&o for A = 0. Upper right, lower left, and lower right: W vs J -1 for A = 0.5, 1, 2, respectively. The 
critical point is where the curves for different system sizes cross. Note that only system sizes much larger than A scale well 
(solid lines), smaller sizes (dashed lines) suffers from finite size corrections. 



To construct a worm one starts at a random lattice 
site Si in space-time. A vortex segment is then created 
in direction a to a new lattice site S2- The direction is 
chosen with the probability 



P a ~ A a = min(l,exp(-AS' (T )), 



(36) 



where AS a is the change in the action S for creat- 
ing a vortex segment in direction a. A new direc- 
tion and lattice site are then chosen from S2- These 
steps are repeated until the worm of vortex segments 
form a closed loop. To ensure detailed balance, the 
final worm is accepted with the probability P = 
min[l, A(noworm)/iV(worm)], where N = A a at the 
initial site, before and after the creation of the loop. Oc- 
casionally the loop constructed this way winds around 
the system, leading to the desired changes in the wind- 
ing number M v . 

The model allows magnetic monopoles in addition to 
the flux lines. The monopoles are updated using a stan- 
dard Metropolis scheme, with trial moves consisting of 
the insertion of a nearest neighbor pair of oppositely 



charged monopoles connected by q flux quanta on the 
link between them. These moves are then accepted with 
probability min(l, exp(— A5)). The creation of a plus 
monopolc on top of a minus monopolc (or vice versa) 
leads to their annihilation so that the procedure de- 
scribed includes creation, destruction, as well as motion 
of monopoles. 



It is straightforward to show that the algorithm de- 
scribed above, combining both the worm updates of the 
flux lines and the Metropolis monopole moves, fulfills 
both ergodicity and detailed balance, and respects the 
constraint V • m = qN. The Wilson loop order parame- 
ter Wn in direction \i is calculated in the simulations as 
(cos(27r M^/q)) where M M is the number of vortex lines 
crossing a plane, perpendicular to direction tt, of the sys- 
tem. 
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FIG. 4: Phase diagram of the model with q = 2. Circles in- 
dicate the points which have been simulated. The parameters 
have been rescaled in order to fit the phase diagram into the 
figure, so that J — > oo occur at the upper side and g — > at 
the right side. The phase diagram connects the 3D XY point 
at g = and the Z2 transition at strong coupling J — > 00. 
The deconfined phase is topologically ordered. 



IV. RESULTS 

For q = 2 we have made simulations for constant 
A = 0, 0.5, 1.0, 2.0, varying J (or g in the A = case) in 
each simulation. The results for W are shown in Fig. [3] 
Below the transition W goes to one with increasing sys- 
tem size. This is the deconfincmcnt phase. Above the 
transition, in the confinement phase, it goes to zero. Pre- 
cisely at the transition W should be independent of sys- 
tem size according to the scaling formula Eq. (|2ip. The 
figures clearly show how W for system sizes much larger 
than A cross at a given J c . However, it is also apparent 
that for large values of A there are strong corrections to 
scaling which influence the smaller sizes. This is not sur- 
prising since the interaction Eq. looks screened only 
at distances much larger than A. It is important to have 
this in mind during the finite size scaling analysis below. 
The phase diagram determined from the simulations is 
shown in Fig. 

As discussed W differs from an ordinary Wilson loop in 
that it does not follow a perimeter law in the deconfined 
phase. In order to demonstrate both area and perimeter 
law, we use a Wilson loop that covers one quarter of the 
system. This configuration is less suitable to locate J c , 
but demonstrates the area and perimeter laws nicely in 
Fig.0 

We also simulate the model with q = 3 to show that 
W works also for larger values of q. The result for W is 
shown in Fig. for A = 0. For q = 3 the model has a first 
order transition when A = 0^ and hence we expect no 
scaling of W. This is consistent with our MC data, since 
we find no good crossing. Instead the transition becomes 



sharper with increasing system size and it is conceivable 
that a discontinuity develops in the thermodynamic limit. 

We now turn to the finite size scaling relation given in 
Eq. (|21|l . and the correlation length exponent v, for q = 2. 
With an appropriate choice of v the data for different 
sizes should collapse onto a single scaling function. This 
is achieved by minimizing 



]T / \w L {x)-W(x) 



dx, 



(37) 



where x = (J 1 — J c 1 )L 1 / 1 ' and W is the mean value 
of Wi,. We calculate v and J c by making a multi pa- 
rameter minimization of Eq. I|37l) . Figure [S] shows two 
of the resulting data collapses. In doing this minimiza- 
tion it is important to make use only of system sizes 
large enough that corrections to scaling are negligible, 
hence we include only those sizes for which well defined 
crossings were found above. We get v = 0.633 ± 0.007, 
0.64 ± 0.015, 0.62 ± 0.03, 0.6 ± 0.15, for A = 0, 0.5, 1, 2, 
respectively. Within error-bars the data for all values of 
A that we have tested are consistent with the expected 
v « 0.63 of the 3D Ising model, which is the dual of 
the Zi gauge theory, i.e., the limiting case obtained for 
A — > 0. This is in sharp contrast to the continuously 
varying exponents found by Sudb0 et alw^ and Smiseth 
et ulM- Their results were obtained from studies of the 
third moment of the action and showed substantial devia- 
tions from the 3D Ising value for A > 0.5. The parameter 
values used in these papers overlap to large extent with 
the present ones, although there is a slight shift in the 
location of the phase diagram resulting from our use of 
the Villain model [Eq. , or equivalently |J3J|] instead of 
the cosine version Eq. . This is not expected to change 
the critical exponents, however. 
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FIG. 5: Area law and perimeter law demonstrated for A — > 0, 
for a Wilson loop W(C) which covers a quarter of a cross 
section of the system and therefore has a finite perimeter and 
area. Dashed and solid lines are exact area and perimeter 
laws, inserted for comparison. 
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FIG. 6: Finite size scaling plot of the W order parameter for A = 0.5 (left) with v = 0.64 and A = 1 (right) with v — 0.62. 
Jc 1 is determined from Fig.^to 1.25 and 2.29, respectively. The inset shows fitting error as a function of v. 



Finally, we may also use W to calculate two quantities 
which are characteristic of topological order, namely the 
gap to vortex excitations E v , and the splitting A of the 
nearly degenerate ground states on a torus. Figure |H1 
shows £E V calculated using Eq. (|34[l as a function of £//?, 
where £ ~ |<5| _i/ . The data nicely follows the scaling 
relation proposed in Eq. (|35|l in the deconfmed phase, 
with a 10. 

The relation between W and A in Eq. (|27|1 allows the 
ground state splitting to be explicitly calculated using 
the Monte Carlo data. An example of this is shown in 
Fig. El where the exponential dependence on system size 
obtained in Eq. (|3L)fl is clearly seen in the deconfinement 
phase (solid line) , while an area law (dashed line) is found 
in the confinement phase. For further examples and dis- 
cussion of this aspect see Ref. 



V. SUMMARY 

In summary, we have studied a compact U(l) gauge 
theory coupled to bosonic matter with gauge charge q us- 
ing duality arguments and Monte Carlo simulations. The 
confinement-deconfinement transition is analyzed using 
a nonlocal order parameter W, which is related to large 
Wilson loops with the loop covering a whole cross section 
of the system, and which directly probes the confining 
properties of the theory. The model can be mapped to 
a dual model with a global q-state clock symmetry, and 
W corresponds to the domain wall energy of that model, 
which suggests a scaling behavior according to Eq. i|21[l . 
We confirm this using Monte Carlo simulations and study 
the details of the deconfinement transition over a range 
of parameters. Our results for the critical exponent v 
support a transition in the 3D Ising universality class for 




0.302 



0.304 



0.306 



0.308 



FIG. 7: W vs g/$l for a model with q = 3 and A = 0. Note 
that this model has a first order transition and no scaling 
behavior is expected. 
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FIG. 8: The energy gap to vortex excitations as a function 
of temperature in the deconfmed phase, for q = 2 and A = 0. 
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FIG. 9: The scaling combination 2A£ 2 /L plotted against 
L/f for q = 2, A = 1. In the topologically ordered decon- 
finement phase the data follows the exponential dependence 
given in Eq. 1301 . as indicated by the solid line. The dashed 
line indicates the area law obeyed in the confinement phase. 



all values of A tested. In particular, we do not see any 
trace of con tinuously varying critical exponents found in 
Refs. I27l28l using a different order parameter. There are, 
however, strong crossover effects in the vicinity of the 
3D XY transition at A — > oo, which make an accurate 
determination of the exponents tricky in this limit. The 
deconfined phase is topologically ordered with gapped 
vortex excitations, with a gap E v ~ £ , which can be 
related to W for a spatial loop. The ground state degen- 
eracy characteristic of topological order is lifted in finite 
sized systems with a gap 2A which is simply related to W 
for a space-time loopi^ This establishes a firm relation- 
ship between topological order and deconfinement and 
allows explicit calculations of the ground state energy 
splitting in the topologically ordered phase of the theory. 
We expect the methods developed here to be useful for 
quantitative studies of topological order in other systems 
as well. 
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APPENDIX: GROUND STATE SPLITTING FOR 
GENERAL g AND q 

In this appendix we give a detailed derivation of the re- 
lation between the ground state splitting and W for mod- 
els with gauge charge q on manifolds with arbitrary genus 
g. The configurations contributing to the partition func- 
tion can be split into separate topological sectors depend- 
ing on the flux ^orrii passing through a space-time cross 

section surrounding hole i. Let Z r r t = (^f e~ f3H f + rri^j 
denote the partition function restricted to the sector la- 
beled by m. Then W % = £~ Z^e 2 ** 1 */*/ £™ Z, Tl . The 
Hamiltonian in the ground state subspace can be written 
as 



f + m)(j 



(A.l) 



where /ig = Eq/2 and h r % = T 7 % (to ^ 0) are tunnel- 
ing amplitudes. Because of translation invariance in the 
index / the eigenstates are 



Jb = 



l 



,2-xik-f/q 



(A.2) 



with corresponding eigenvalues 

E% = h^2 cos ^27rfc ■ rh/c 



(A.3) 



To relate the energies to the W? that are easily calculated 
in a simulation we use 



^ = (kle-W \k) = 



(A.4) 



q 2g /Z(f 

If' 



f)e 



27Tik-(f-f')/q 



Z^e 2 ™*-™'* = W % Za, 

rh m 

from which Eq. (|3 1|) follows. 
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